# load libraries ----
# reporting
library(knitr)
# visualization
library(ggplot2)
library(ggthemes)
library(patchwork)
library(ggrepel)
# data wrangling
library(dplyr)
library(tidyr)
library(readr)
library(skimr)
library(janitor)
library(magrittr)
library(lubridate)
# lake models
library(GLMr)
library(glmtools)
# mapping
library(sf)
library(ggspatial)
library(rnaturalearth)
library(rnaturalearthhires)
library(rnaturalearthdata)
# set other options ----
options(scipen=999)
knitr::opts_chunk$set(
tidy = FALSE,
message = FALSE,
warning = FALSE)31 Teleconnections I
Learning Objectives
After completing this tutorial you should be able to
- describe the concepts of macrosystems and teleconnections and how different ecological processes can interact at local, regional, and global scales.
- set up and run lake models to simulate lake temperatures, thermal structure, and ice cover in multiple lakes.
- formulate hypotheses about the effects of teleconnected climate scenarios on different lake modules and test them using simulations.
- describe how local characteristics modify global-scale climate forcing effects on lake temperatures and ice cover.
Download the directory for this project here, make sure the directory is unzipped and move it to your bi328 directory. You can open the Rproj for this module either by double clicking on it which will launch Rstudio or by opening Rstudio and then using File > Open Project or by clicking on the Rproject icon in the top right of your program window and selecting Open Project.
There should be a file named 31_elnino_I.qmd in that project directory. Use that file to work through this tutorial - you will hand in your rendered (“knitted”) quarto file as your homework assignment. So, first thing in the YAML header, change the author to your name. You will use this quarto document to record your answers. Remember to use comments to annotate your code; at minimum you should have one comment per code set1 you may of course add as many comments as you need to be able to recall what you did. Similarly, take notes in the document as we discuss discussion/reflection questions but make sure that you go back and clean them up for “public consumption”.
1 You should do this whether you are adding code yourself or using code from our manual, even if it isn’t commented in the manual… especially when the code is already included for you, add comments to describe how the function works/what it does as we introduce it during the participatory coding session so you can refer back to it.
Before you get started, let’s make sure to read in all the packages we will need.
31.1 Teleconnections
We have already seen how local and regional drivers can impact the thermal structure of lakes. Local impacts could include something like land cover change which alters the stream inflow/outflow of a lake impacting temperature. At a regional scale climate change can alter weather patterns, including both extreme events which disrupt regular seasonal patterns or also gradual change which raises lake temperatures overall.
Next, we are going to explore how ecosystems can be influenced on a global scale by teleconnections, i.e. meteorological, societal, and/or ecological phenomenon that link remote regions via cause and effect relationships. Simulations are especially powerful in this context, because it allows researchers to manipulate drivers independently from each other to disentangle how they impact ecosystems. By running different teleconnection scenarios, scientists are then able to better predict how ecosystems will respond to remote drivers.
In our module, we will specifically look at lake thermal structure and ice cover and determine how it is impacted by climate teleconnections to decadal events such as the El Niño/Southern Oscillation (ENSO). This is an example of a global driver. ENSO period with warmer ocean surface temperatures in the Pacific Ocean affect regular trade wind patterns. This altered atmospheric circulating impacts air temperature and precipitation patterns globally. How ENSO ends up impacting a particular lake will be further mediated by local and regional dynamics.
Our central question is:
How do global ENSO teleconnections interact with local & regional drivers to affect lake temperatures, thermal structure, and ice cover?
We can address this question using the General Lake Model we used to understand the effect of local & regional drivers on the thermal structure of lakes. To explore those interactions, we varies the local/regional drivers (weather/climate) and ran different scenarios (described in our met file) but kept the lake we were using constant, i.e. the morphology of the lake and other specific characteristics described in the nml file did not change.
ENSO has large-scale effects, however, the way an El Nino year plays out will vary by region, additional local effect will further modify the impacts on the thermal structure of each lake.
Consider this
During El Nino years there is a shift in the climate pattern of the Pacific Ocean, resulting in unusually warm surface waters in the the eastern equatorial Pacific Ocean. Use the figure above to describe regional differences in climate patterns in North America are impacted.
To be able to explore our question we will need nml files describing lakes in different regions across North America. You should a folder called sim_lakes to your project directory. If you take a look in there, you fill find individual folders for a set of lakes, in each folder you have an nml file describing that lake along with a met file containing observed meteorological data. These are all lakes that are part of the GLEON network.
Let’s take a look at where those lakes are situated.
lakes <- read_delim("data/lake_characteristics.txt", delim = "\t") %>%
clean_names()
kable(lakes)| lake_name | state | latitude | longitude | elevation_m | max_depth_m | surface_area_km2 | trophic_state | mixing_regime | watershed_land_use | neon_domain | data_provider | links |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Barco Lake | Florida | 29.67594 | -82.00823 | 28 | 7 | 0.12 | Oligotrophic | Warm monomictic | Forest | 3 | NEON Core Aquatic | http://www.neonscience.org/field-sites/field-sites-map/BARC |
| Crampton Lake | Wisconsin | 46.21025 | -89.47345 | 591 | 18 | 0.26 | Oligotrophic | Dimictic | Forest | 5 | NEON Core Aquatic | http://www.neonscience.org/field-sites/field-sites-map/CRAM |
| Falling Creek Reservoir | Virginia | 37.30333 | -79.83722 | 507 | 9 | 0.12 | Eutrophic | Dimictic | Forest | 7 | GLEON; Carey Lab | https://www.westernvawater.org/drinking-water/water-sources-and-treatment |
| Lake Mendota | Wisconsin | 43.10970 | -89.42060 | 259 | 25 | 39.39 | Eutrophic | Dimictic | Agricultural/urban | 6 | LTER- North Temperate Lakes | https://lter.limnology.wisc.edu/researchsite/lake-mendota |
| Lake Sunapee | New Hampshire | 43.38020 | -72.05270 | 333 | 33 | 16.74 | Oligotrophic | Dimictic | Forest | 1 | GLEON; Lake Sunapee Protective Association | http://www.lakesunapee.org/lake-sunapee-history/ |
| Prairie Pothole | North Dakota | 47.13012 | -99.25302 | 565 | 3 | 0.11 | Mesotrophic | Dimictic | Grasslands | 9 | NEON Core Aquatic | http://www.neonscience.org/field-sites/field-sites-map/PRPO |
| Suggs Lake | Florida | 29.68770 | -82.01785 | 29 | 6 | 0.73 | Mesotrophic | Warm monomictic | Forest | 3 | NEON Core Aquatic | http://www.neonscience.org/field-sites/field-sites-map/SUGG |
| Toolik Lake | Alaska | 68.62956 | -149.61051 | 738 | 26 | 1.46 | Oligotrophic | Dimictic | Tundra | 18 | NEON Relocatable Aquatic; LTER- Arctic | http://arc-lter.ecosystems.mbl.edu/toolik-main |
Table 31.1: Characteristics for eight lakes in the GLEON network.
Let’s make a map displaying their geographic locations. We’ve previously used geom_text() to add text (labels) to maps. But those weren’t exactly pretty maps. We are going to add a new tool to our toolbox for generating maps, geom_label_repel(). This function works especially well if you want to combine using a shape to plot a point and add a label that is placed in a way that it does not overlap with that point or other labels on the plot.
# dowhload world map as shapefile
world <- ne_countries(scale = "medium", returnclass = "sf")
# download us states as suhape file
us_states <- ne_states(country = "United States of America", returnclass = "sf")
# lat/long for map extent
x_min <- min(lakes$longitude) - 2
x_max <- max(lakes$longitude) + 2
y_min <- min(lakes$latitude) - 2
y_max <- max(lakes$latitude) + 2
# set color for fill
map_color <- "khaki3"
# create plot
ggplot() +
geom_sf(data = world, color = "black", fill = map_color) + # plot outline of countries
geom_sf(data = us_states, color = "black", fill = NA) + # plot outline of states
coord_sf(xlim = c(x_min, x_max),
ylim = c(y_min, y_max)) + # plot boundaries for map
geom_point(data = lakes, aes(x = longitude, y = latitude), # add sites
size = 2) +
geom_label_repel(data = lakes,
aes(x = longitude, y = latitude, label = lake_name),
size = 3)
Don’t worry - you won’t have to run simulations for every lake - we will divide and conquer!
31.2 Set up simulation
If you take a look at the sim_lakes folder, you should see that it contains individual folders for a series of lakes. In this example, we will use Lake Sunapee - this means that when you run the code for your lake, you will need to be sure to change the lake name accordingly.
First, let’s create a vector with our lake name2.
2 The advantage for doing it this way instead of directly in the individual code lines is that we only have to set it once, this makes the code a lot more reusable because we don’t have to dig through the code and find all the arguments that set a file path or specify the lake name.
LakeName <- "Sunapee"Next, we will set the file path for the sim folder. Each lake sim folder contains the input files needed to run GLM: the nml file which describes the lake, the met file describing the drivers of thermal structure, and input files describing inflow and outflow to the lake.
sim_folder <- file.path("sim_lakes", LakeName)Next, we’ll specify the file path for the nml file, read it in, and take a look at it.
# define file path
nml_file <- file.path(sim_folder, "glm2.nml")
# read nml file
nml <- read_nml(nml_file)
# view file
print(nml)&glm_setup
sim_name = 'MacrosystemsEDDIE_Sunapee'
max_layers = 250
min_layer_vol = 0.025
min_layer_thick = 0.05
max_layer_thick = 0.75
Kw = 0.17
coef_mix_conv = 0.08213
coef_wind_stir = 1.27
coef_mix_shear = 0.296
coef_mix_turb = 0.526
coef_mix_KH = 0.415
coef_mix_hyp = 0.453
/
&morphometry
lake_name = 'Sunapee'
latitude = 43.395
longitude = -72.053
bsn_len = 12870
bsn_wid = 3000
bsn_vals = 82
H = 299.43, 299.943, 300.443, 300.943, 301.443, 301.943, 302.443, 302.943, 303.443, 303.943, 304.443, 304.943, 305.443, 305.943, 306.443, 306.943, 307.443, 307.943, 308.443, 308.943, 309.443, 309.943, 310.443, 310.943, 311.443, 311.943, 312.443, 312.943, 313.443, 313.943, 314.443, 314.943, 315.443, 315.943, 316.443, 316.943, 317.443, 317.943, 318.443, 318.943, 319.443, 319.943, 320.443, 320.943, 321.443, 321.943, 322.443, 322.943, 323.443, 323.943, 324.443, 324.943, 325.443, 325.943, 326.443, 326.943, 327.443, 327.943, 328.443, 328.943, 329.443, 329.943, 330.443, 330.943, 331.443, 331.943, 332.343, 332.443, 332.543, 332.643, 332.743, 332.843, 332.943, 333.043, 333.143, 333.243, 333.343, 333.443, 333.543, 333.643, 333.743, 333.943
A = 1, 16.90309, 38.87712, 64.23176, 158.88908, 7254.80801, 25205.89398, 33207.81875, 41551.18602, 51853.6219, 62734.1436, 79030.41666, 105037.51731, 133941.80829, 168576.24818, 214836.63623, 269137.82616, 327112.05845, 401221.9844, 483528.22071, 575947.57823, 676930.04325, 797467.69789, 925704.7119, 1088626.8751, 1283336.99713, 1505854.39934, 1740820.93038, 1988527.32348, 2234911.89462, 2485285.59597, 2723309.89675, 2978076.71227, 3264022.97524, 3653258.97559, 4050485.06854, 4414109.50064, 4796077.1705, 5158332.45198, 5579849.98152, 6016348.72339, 6468883.43068, 6928267.27171, 7476296.00935, 8028820.97004, 8623779.45825, 9274820.72252, 9852960.63206, 10365149.7392, 10843023.87491, 11253594.96068, 11644117.28649, 12013032.38706, 12356310.56475, 12673823.35603, 12988489.66626, 13285122.06535, 13576406.32546, 13863498.61821, 14148319.13512, 14418854.84716, 14691317.51194, 14951571.07182, 15234317.57908, 15510517.51797, 15780918.00526, 16009460, 16071489, 16135107, 16203453, 16273655, 16349525, 16489737, 16603406, 16760708, 16788809, 16804233, 16826128, 16844170, 16863779, 16885473, 16934251.6
/
&time
timefmt = 2
start = '2012-11-01 00:00:00'
stop = '2013-12-31 23:00:00'
dt = 3600
timezone = -5
/
&output
out_dir = '.'
out_fn = 'output'
nsave = 24
/
&init_profiles
lake_depth = 33.7
num_depths = 6
the_depths = 0, 4, 8, 12, 16, 20
the_temps = 11, 11, 11, 11, 11, 11
the_sals = 0, 0, 0, 0, 0, 0
/
&meteorology
met_sw = .true.
lw_type = 'LW_IN'
rain_sw = .false.
atm_stab = .false.
catchrain = .false.
rad_mode = 2
albedo_mode = 1
cloud_mode = 4
meteo_fl = 'met_hourly.csv'
subdaily = .true.
wind_factor = 1
sw_factor = 1
lw_factor = 1
cd = 0.0013
ce = 0.0013
ch = 0.0013
rain_threshold = 0.01
runoff_coef = 0.3
/
&inflow
num_inflows = 1
names_of_strms = 'stream'
strm_hf_angle = 65
strmbd_slope = 3
strmbd_drag = 0.016
inflow_factor = 1
inflow_fl = 'inflow.csv'
inflow_varnum = 3
inflow_vars = 'FLOW','SALT','TEMP'
/
&outflow
num_outlet = 1
flt_off_sw = .false.
outl_elvs = 331
outflow_fl = 'outflow.csv'
outflow_factor = 1.05
seepage_rate = 0
/
&snowice
snow_albedo_factor = 1
snow_rho_max = 500
snow_rho_min = 100
/
&sed_heat
sed_temp_mean = 12
sed_temp_amplitude = 11
sed_temp_peak_doy = 242.5
/One of the arguments in the nml file specifies the file with the information on drivers (i.e. it directs to the met_hourly.csv) file. Let’s take a look at the meteorological input data for the duration of the simulation.
plot_meteo(nml_file) 
31.3 Run & analyze baseline simulation
Now that we have everything set up, we can run the model.
run_glm(sim_folder, verbose = TRUE)[1] 0To analyze the output, let’s specify the file path for the output file that is now in our lake simulation folder (output.nc). Again, this will be our baseline data set based on observed meteorological data for a specific year (in our example 2013).
baseline <- file.path(sim_folder, "output.nc")Let’s take a look at the variables that were output as part of our GLM simulation run.
# pull variable names from output
var_names <- sim_vars(baseline)
# print variables in output
kable(var_names,
caption = "Variable names and units for lake simulation.") | name | unit | longname |
|---|---|---|
| evap | m/s | evaporation |
| extc_coef | unknown | extc_coef |
| hice | meters | Height of Ice |
| hsnow | meters | Height of Snow |
| hwice | meters | Height of WhiteIce |
| I_0 | 10E-6m | Shortwave |
| NS | Number of Layers | |
| precip | m/s | precipitation |
| rad | unknown | solar radiation |
| rho | unknown | density |
| salt | g/kg | salinity |
| temp | celsius | temperature |
| Tot_V | m3 | lake volume |
| u_mean | m/s | mean velocity |
| u_orb | m/s | orbital velocity |
| V | m3 | layer volume |
| wind | m/s | wind |
Variable names and units for lake simulation.
One of the variables we are especially interested in is hice which describes the lake thickness on the surface of the lake.
Part of our central question is understanding how water temperatures at the lake surface and bottom (below the thermocline) be affected by El Nino teleconnections.
To estimate the bottom depth we can extract the lake depth in meters from the output.
# extract lake depth
LakeDepth <- get_surface_height(baseline)
head(LakeDepth) DateTime surface_height
1 2012-11-02 33.70000
2 2012-11-03 33.69793
3 2012-11-04 33.69418
4 2012-11-05 33.69017
5 2012-11-06 33.68645
6 2012-11-07 33.68364We are going to want to extract the lake temperature data from our baseline simulation into a data.frame so we can compare it to our teleconnection scenarios down the line. We will extract it at the surface and for the bottom by setting the argument z_out to 0 (surface) and the minimum observed depth over the simulation (min(LakeDepth$surface_height)).
The function automatically names the temperatures after the depth - so you will want to rename them something sensible. You will need to pull up the original column names for your lake, and then you will need to change the column names in the rename() function accordingly!
# pull temperature data
LakeTemp <- get_temp(file = baseline, reference = "surface",
z_out = c(0, min(LakeDepth$surface_height)))
# get column names
colnames(LakeTemp)[1] "DateTime" "temp_0" "temp_33.3833669399235"# rename columns
LakeTemp <- LakeTemp %>%
rename(Baseline_SurfaceTemp = temp_0,
Baseline_BottomTemp = temp_33.3833669399235)
Pointers
Hint - to plot both temperatures in the same plot, you will need to first create a tidy data set!
If you look at the figure below, you can see that you will need a column for date (x-axis), temperature (y-axis) and and column that contains simulated scenarios (this becomes the color coding for your line graphs).
Once you do that, it should be straightforward to create a plot that look something like this:
Let’s plot a heatmap of our lake thermal structure for the baseline scenario.
plot_temp(file = baseline, fig_path = FALSE) 
Next, we will extract the ice cover on the lake for our baseline simulation.
LakeIce <- get_var(baseline, "hice") %>%
rename(ice_baseline = hice)31.4 Give it a whirl
Visualize the change of surface ice for your lake over the course of the year and then describe your results. Remember to be specific - when is there ice? How thick is it? Include a title & label your axes!
:::
Solution
Your plot should look something like this:
Now that we have a pretty good idea what our baseline data set looks like in terms of the thermal structure in a normal year, our next step will be to simulate what the thermal structure will look like in a typical or extreme El Nino year compared to an extreme El Nino year.
31.5 Acknowledgments
These activities are based on the EDDIE Teleconnections module.3
3 Farrell, K.J. and C.C. Carey. 18 May 2018. Macrosystems EDDIE: Teleconnections. Macrosystems EDDIE Module 3, Version 1.